1. Field of the Invention
This invention relates to a method for interpolating between measured waveforms obtained from an acoustic, elastic or other experiment which obeys the wave equation. Interpolation between measured waveforms is accomplished by successive forward and backward applications of wave equation datuming. This invention further relates to seismic exploration and to a method for seismic trace interpolation by applying wave equation datuming techniques.
2. Description of the Prior Art
Physical experiments are often conducted in which recorded measurements which pertain to the outcome of the physical experiment obey the wave equation. Most often, the measurements recorded in such an experiment are related to the time history of the particle motion, displacement, acceleration or pressure at the location of a transducer. For example, in the field of underwater acoustics, a source may be used to generate sound waves into the ocean. An array of hydrophones towed behind a boat will be used to record the sound waves impinging upon the array. The recorded sound waves, which may be used to provide information regarding the location of objects beneath the ocean surface, will obey the wave equation. In the field of exploration seismology, an example of a physical experiment which obeys the wave equation would be the generation of seismic energy by a source and the use of an array of geophones spread out on the surface of the earth for recording elastic waves which are reflected from the subsurface and which impinge upon the array.
In seismic exploration, it is common practice to deploy a large array of geophones on the surface of the earth and to record the vibrations of the earth at each geophone location to obtain a collection of seismic traces. The traces are sampled and recorded for further processing. When the vibrations so recorded are caused by a seismic source activated at a known time and location, the recorded data can be processed by a computer in known ways to produce an image of the subsurface. The image thus produced is commonly interpreted by geophysicists to detect the possible presence of valuable hydrocarbons.
Seismograms are commonly recorded as digital samples which represent the amplitude of a received seismic signal as a function of time. Since seismograms are usually obtained along a line of exploration on the surface of the earth, the digital samples can be formed into x-t arrays with each sample in the array representing the amplitude of the seismic signal as a function of horizontal distance and time. When such arrays are visually reproduced, by plotting or the like, seismic sections are produced. A seismic section depicts the subsurface layering of a section of the earth. It is the principal tool which the geophysicist studies to determine the nature of the earth's subsurface formations. Before an array of seismic samples or traces can be converted into a seismic section for interpretation by geophysicists, the array must be extensively processed to remove noise and to make reflection events discernable.
It frequently occurs, however, that insufficient data is recorded data for proper interpretation of the characteristics of the subsurface formation. The most common data insufficiency problem is that the acquired data includes either seismic traces with no recorded data or seismic traces which clearly contain severe noise contamination such as spikes. These and other insufficiencies in the acquired data may be traced to any one of numerous causes. For example, a data insufficiency may result when seismic traces are recorded in the wrong location due to intrabed reflections of the generated seismic energy. Another typical data insufficiency arises when it is determined after seismic exploration has been completed that the collected traces need to be interpolated to a finer spatial sampling for proper interpretation. Here, the selection of too large a spacing between geophones prior to commencing exploration is typically the cause. A seismic trace without data may be produced when one or more geophones fail during exploration. Seismic traces severely contaminated with noise may be inadvertently produced due to multiple reflections or ground roll.
Too often, cost or time restraints will prevent the reshooting of the unsatisfactory data. Under such circumstances, the most acceptable option available to the geophysicist who must complete an evaluation of the formation based upon the insufficient data, is to estimate the missing or bad traces from the set of collected traces.
There are numerous methods for determining estimates of data points from measured physical data which may be classified as interpolation methods. For example, many methods exist for the interpolation of functions which are dependent on only one parameter such as time. These include both the well known linear interpolation methods and the equally well known finite difference interpolation methods. Linear interpolation is based upon the proposition that for an argument "x", the value of a dependent variable "y" will vary linearly over the interval under investigation. Finite difference interpolation methods are based upon the proposition that given "y=f(x)", inverse interpolation may be used in a process for finding "x" at a value of "y" intermediate to two tabulated values. For finite difference interpolation methods, any one of several alternate procedures such as Lagrange's formula, successive approximation, or reversion of series may be used to obtain the desired result.
Standard practice among geophysicists faced with seismic traces which contain either no recorded data or which are severely contaminated has been to exclude such traces, commonly referred to as "null" traces, from the otherwise satisfactory data set. The collected seismic data would be processed using conventional methods without the excluded data. When the missing trace was necessary for proper processing of the seismic data, prior attempts to restore the missing trace and create a trace with events consistent with nearby coherent events focused upon combining traces near the missing trace in the x-t domain to create the missing trace. The use of such prior methods restricted conventional multi-dimensional trace interpolation upon the apparent characteristics of the measured traces and did not permit trace interpolation to use the characteristics of the underlying physical experiment being performed.
The above mentioned trace interpolation methods are inherently 1-dimensional, for example, f(x)=b. These methods may not be applied to data sets which are 2 or 3-dimensional, such as f(x,t)=A, or f(x,y,t)=B. For example, a seismic trace measures amplitude as a function of time. Thus, each seismic trace in the gather is used to analyze subsurface formations which are separated by a distance dependent on the location of the geophones along the surface. This spatial variation may be only in one, two or three directions, depending on the location of the geophones utilized in the gather.
Alternate methods to the interpolation methods described above were developed which address the problems of interpolating functions having spatial and temporal dependencies. Too often, however, such alternate methods are based upon the observed statistics of the trace gather rather than upon the physics from which the experiment came. One such conventional solution which relies upon the statistics of the trace gather and is popular in the field of seismic exploration utilizes t frequency-wavenumber transform approach.
In the processing of seismograms, x-t arrays have been transformed into arrays representing amplitude as a function of frequency and wave number. This is commonly referred to as a "frequency-wave number" or "f-k" transformation. The f-k transformation has been used as a tool to study seismic filtering. F-k transforms are routinely used to represent data collected by a large array of sensors. Usually, the f-k representations are computed by Fast Fourier Transforms (hereafter referred to as FFTs) which transform the data from the time and spatial domain to the frequency and wavenumber domain by utilizing the corresponding multi-dimensional Fourier transform. The resulting data representations are parameterized by frequencies, wave numbers (spatial frequencies), amplitudes, and phases. More specifically, for each frequency there is a collection of wave numbers, and for each frequency-wave number pair there is a complex number representative of an amplitude and a phase. The data is then inverse Fast Fourier transformed back to the time and spatial domain. After this inverse transform, the traces have a different spatial sampling and possible a different temporal sampling.
In U.S. Pat. No. 4,594,693 issued to Pann et al, seismic trace interpolation is carried out by inserting zero amplitude traces between the seismic traces in a section where spatial aliasing is a problem. The traces are then transformed into an f-k array. The f-k array is filtered with a filter which rejects samples in a region of frequency and wave number which exhibits aliasing. The filtered f-k array is then transformed back into a seismic section representing amplitude as a function of time and distance. One limitation on the use of such methods, however, it that seismic interpolation using f-k transforms requires that the measured traces be uniformly sampled in space. This requirement of a regular trace spacing limits the usefulness of f-k transforms for seismic interpolation applications.
Finally, it is known to use a technique generally referred to as "velocity replacement" to remove the distortion of information gathered during seismic exploration which is caused by changes in the velocity of seismic energy between irregular, i.e. non-planar, layers in the subsurface formation being explored. Due to the varied characteristics of the layers, seismic energy travels at different velocities in each layer. To remove distortion from seismic data recorded at the surface which has traveled through such an irregular subsurface layer, the seismic data is propagated at a first velocity to a first real datum, corresponding to the bottom of the irregular layer, and then propagated at a second velocity to a second real datum, corresponding to the top of the layer. Further, as each layer in the surface has different characteristics which affect the velocity of seismic energy passing through that layer, the actual velocity of the irregular layer is used as the first velocity, and the velocity of the medium below the layer is used as the second velocity. For this reason, velocity replacement methods have consistently used different velocities for each propagation.